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I.    INTRODUCTION 
A.   THE  IDEA  OF  ARMA  MODELING 

The  goal  of  linear  modeling  is  to  accurately  represent  an  observed  data  sequence 
as  the  output  of  a  linear  filter.  The  idea  of  representing  a  complicated  process  with  a 
comparatively  simpler  model  has  many  different  applications.  Curve  fitting  in  math- 
ematical modeling,  analysis  of  electronic  devices  using  equivalent  circuits,  and  system 
transfer  functions  in  automatic  control  are  just  a  few  examples  outside  the  area  of  dig- 
ital signal  processing.  Parametric  modeling  has  also  a  large  number  of  applications  in 
signal  processing.  Currently  there  is  considerable  interest  in  the  parametric  modeling 
approach  to  spectral  estimation.  In  speech  processing  the  applications  include  digital 
transmission,  storage,  and  synthesis  of  the  speech  signal.  Our  particular  interest  is 
the  modeling  of  sonar  signals,  such  as  biologies  and  other  underwater  acoustic  data. 
This  work  forms  part  of  an  overall  research  program  in  sonar  signal  modeling.  The 
research  will  help  to  understand  the  relative  benefits  of  signal  domain  algorithms 
versus  algorithms  based  on  coefficients  of  the  transfer  function.  It  is  hoped  that  the 
method  developed  in  this  thesis  will  become  an  important  tool  in  the  overall  effort 
for  sonar  signal  modeling. 

In  linear  modeling  the  filter  used  to  generate  the  data  sequence  is  usually  rep- 
resented by  a  linear  difference  equation  with  constant  coefficients.  The  2-transform 
of  this  type  of  system  is  a  rational  polynomial  function.  Three  type  of  models  are 
derived  from  this  kind  of  systems;  they  are  known  as  autoregresive  (AR),  moving 
average  (MA),  and  autoregresive  moving  average  (ARMA). 

Much  work  has  been  done  on  AR  models,  which  correspond  to  all-pole  systems. 
The  reason  for  that  is  that  the  parameters  for  the  model  can  be  obtained  by  solving 
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linear  equations,  and  a  great  body  of  theory  has  been  developed  that  applies  to 
this  problem  [Ref.  1,  2].  Relatively  much  less  work  has  been  done  with  MA  and 
ARMA  models.  However,  since  MA  models  have  limited  applications  and  are  almost 
as  difficult  to  obtain  as  ARMA  models,  most  interest  centers  on  the  latter.  The 
filter  in  this  type  of  models  has  both  poles  and  zeros,  and  the  design  fundamentally 
involves  nonlinear  equations.  A  properly  designed  ARMA  model  can  provide  better 
performance  than  an  AR  model,  with  a  smaller  number  of  parameters.  ARMA 
modeling  is  the  topic  of  this  thesis. 

B.    WHY  AN  ITERATIVE  ALGORITHM  IN  THE  SIG- 
NAL DOMAIN 

There  are  many  different  approaches  to  the  problem  of  ARMA  modeling.  The 
majority  of  them  are  based  on  statistical  techniques  [Ref.  3].  Some  of  these  methods 
regard  the  data  as  a  realization  of  a  random  process  while  others  focus  on  the  data  as 
given  [Ref.  1].  Data  oriented  methods  try  to  minimize  some  criterion  that  estimates 
how  well  the  model  fits  the  data,  in  most  cases  the  least  squares  error  between  the 
signal  and  the  model.  Stochastic  approaches  may  attempt  to  estimate  the  model 
parameters  directly  from  the  data  by  solving  nonlinear  equations  or  by  spectral  fac- 
torization. The  maximum  likelihood  procedure,  for  example  [Ref.  1,  4],  is  essentially 
nonlinear.  A  number  of  indirect  methods  have  been  developed  that  modify  the  norm 
of  the  error  by  separating  the  AR  and  MA  parts  of  the  problem  so  that  at  least  some 
of  the  equations  to  estimate  the  parameters  become  linear.  This  approach  is  found 
in  procedures  such  as  the  Prony's  method,  Shank's  method,  and  the  least  squares 
modified  Yule- Walker  method  [Ref.  1,  2,  5,  6].  A  different  type  of  approach  replaces 
the  nonlinear  problem  with  iteration  while  trying  to  solve  for  the  AR  and  MA  pa- 
rameters simultaneously.  The  iterative  prefiltering  method  of  Steiglitz  and  McBride 
[Ref.  7,  8]  is  of  this  type. 


Our  method  is  also  of  the  latter  type.  However,  the  advantage  is  that  it  works 
directly  with  the  poles  of  the  rational  model,  which  we  know  affect  the  performance 
of  the  system.  The  poles  are  displaced  in  specific  directions  so  that  the  new  model 
minimizes  the  error  between  the  model  output  and  the  original  signal.  Iterative 
prefiltering  on  the  other  hand  works  with  the  coefficients  of  the  transfer  function, 
so  it  is  difficult  or  impossible  to  predict  its  effects  on  the  poles  and  zeros  of  the 
system.  The  new  algorithm  is  much  more  dependable  with  respect  to  convergence 
than  the  iterative  prefiltering  algorithm  because  it  moves  poles  and  zeros  specifically 
to  minimize  the  error  between  the  model  and  the  original  signal. 

C.   THESIS  OUTLINE 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  II  presents  the 
modeling  methods  that  are  used  in  this  thesis  and  gives  a  brief  explanation  of  all 
of  them.  Chapter  III  introduces  the  reader  to  the  theory  of  multidimensional  op- 
timization by  gradient  methods  and  develops  the  iterative  Prony  method,  which  is 
the  main  contribution  of  this  thesis.  Chapter  IV  presents  the  results  of  testing  the 
algorithm  on  simulated  and  real  acoustic  data  and  compares  these  results  with  those 
obtained  using  iterative  prefiltering.  Chapter  V  gives  conclusions  and  suggestions  for 
future  research. 


II.    ARMA  MODELING  OF  SIGNALS 

A.   MODELING  METHODS  USED  IN  THIS  THESIS 

This  thesis  deals  with  deterministic  approaches  to  ARMA  modeling.  Two  types 
of  modeling  methods  are  considered.  First  we  have  non-iterative  methods  like  Prony's 
method  and  its  alternate  signal  domain  form  [Ref.  1];  second  we  consider  iterative 
methods  like  iterative  prefiltering  and  the  new  iterative  Prony  method  developed  in 
this  thesis. 

The  goal  of  Prony's  method  is  to  represent  a  given  sequence  x[n]  as  the  impulse 
response  of  a  linear  time  invariant  (LTI)  system.  In  the  transform  (z)  domain  this 
representation  has  the  form 

XW-lg  (2.1) 

where  X(z)  is  the  r-transform  of  x[n]  and  B(z)/A(z)  represents  the  transfer  function 
of  the  system.  This  approximation  is  explained  in  more  detail  in  the  next  section. 
What  has  become  known  as  Prony's  method  in  the  current  signal  processing  literature 
differs  from  Prony's  original  work  in  some  respects.  Our  basic  form  of  Prony's  method 
solves  for  the  coefficients  of  the  transfer  function  in  (2.1). 

The  signal  domain  form  of  Prony's  method  is  closer  to  Prony's  original  work 
[Ref.  1,  p. 560]  and  seeks  to  represent  the  data  in  terms  of  a  set  of  damped  exponen- 
tials as 

x[n]  «  cxr1  +  c2r™  + \-  cPrP 

where  the  r*  are  the  roots  of  the  denominator  polynomial  A(z)  and  the  Ck  are  the 
complex  coefficients  required  for  the  expansion.  Both  forms  involve  linear  equations 
and  least  squares  techniques. 
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Figure  2.1:  Block  diagram  for  the  direct  method  for  signal  modeling. 

The  iterative  prefiltering  method  used  is  due  to  Steiglitz  and  McBride  [Ref.  7,  8] 
and  attempts  to  match  the  data  with  a  model  of  finite  order  using  an  efficient  iterative 
approach.  It  differs  from  Prony's  method  in  that  it  solves  for  both  numerator  and 
denominator  polynomial  coefficients  simultaneously  at  each  iteration. 

Whenever  a  signal  is  modeled  using  a  fixed-order  rational  polynomial  model,  an 
approximation  has  to  be  made  and  some  kind  of  measure  has  to  be  used  to  determine 
the  "goodness"  of  the  model.  In  this  thesis  the  least  squares  error  norm  (rn6ox/2 
norm)  is  used  to  measure  the  approximation  error.  This  norm  measures  the  energy 
of  the  error  and  is  the  norm  most  widely  used  primarily  due  to  its  mathematical 
tractability  [Ref.  2]. 

B.   OVERVIEW  OF  MODELING  METHODS 
1.    Prony's  Method 

The  derivations  of  all  the  modeling  methods  presented  in  this  section  follow 
those  in  [Ref.  1,  pp.  550-564]  and  [Ref.  2].  As  stated  above,  Prony's  method  aims 
at  representing  a  signal  as  the  impulse  response  of  an  LTI  system.  Figure  2.1  shows 
an  implementation  of  this  approximation  which  is  known  as  the  direct  method.  The 
system  function  X(z)  in  the  transform  domain  is  a  rational  polynomial  function 
B(z)/A(z)  with  Q  zeros  and  P  poles.  The  error  signal  e[n]  is  computed  as  the 
difference  between  the  response  of  the  system  x[n]  and  the  given  signal  x[n],  i.e. 

e[n]  =  x[n]  -  x[n].  (2.2) 


The  LTI  system  is  chosen  to  minimize  the  sum  of  squared  errors 


$D  =  I^|e[™]|2- 


(2.3) 


This  problem  leads  to  nonlinear  equations  whose  solution  (if  a  unique  solution  actu- 
ally exists)  turns  out  to  be  a  very  difficult  task  [Ref.  1.2].  To  avoid  these  difficulties, 
a  number  of  indirect  methods  have  been  developed  for  modeling.  Prony's  method  is 
one  of  these  procedures  and  can  be  derived  as  follows.  The  LTI  system  satisfies  the 
difference  equation 


x[n]  +  a\x{n  —!]  +  •••  +  apx[n  —  P]  —  b06[n]  +  &1<5[n  —  !]  +  ■••  +  &Q<5[rc  —  Q]. 


(2.4) 


If  the  requirement  that 


r[n]  =  x[n],         n  =  0, 1, . . .  ,NS  —  1, 


is  applied  to  (2.4),  where  Ns  is  the  length  of  the  data,  and  the  difference  equation  is 
evaluated  for  n  =  0, 1, . . .  ,Na  —  1,  the  result  is  the  matrix  equation 
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where 


XR  = 


x[0] 
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x[l] 

x[Q] 

0 

x[2] 

x[l] 

x[0] 

an( 


X4  = 


x[Q]    x[Q-\]    x[Q-2] 
x[Q  +  l]      x[Q]  x[Q-l] 

x[Ns-l]    x[Ns-2]    x[Ns-S] 


x[Q  -  P] 
x[Q-P+l] 

x[Ns  -  P  -  1] 


(2.7) 


(2.8) 


and  a  and  b  are  the  vectors  of  coefficients  appearing  in  (2.5).  The  lower  partition 


X4a  =  0 


(2.9; 


represents  an  overdetermined  set  of  linear  equations  that  need  not  have  an  exact 
solution.  This  set  of  equations  can  be  solved  by  least  squares,  where  a  is  chosen  to 
minimize  the  least  squares  norm  of  the  equation  error  Sa  =  ||eA||2  in 


X4a  =  eA. 


(2.10) 


This  leads  to  a  set  of  linear  equations  called  the  normal  equations,  which  can  be 
written  compactly  as  [Ref.  1,  pp. 536-537] 


X*ATXA )  a  = 


sA 

0 


(2.11) 


and  can  be  solved  for  a  and  Sa-  Once  the  vector  a  is  known,  it  can  be  substituted 
in  the  upper  partition  of  (2.6) 

b  =  XBa  (2.12) 

to  solve  for  the  vector  b.  Although  it  is  referred  here  simply  as  Prony's  method, 
this  procedure  is  also  known  as  the  modern  Prony  method  or  the  extended  Prony 
method. 
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Figure  2.2:  Block  diagram  for  the  indirect  modeling  problem. 

Although  Prony's  method  is  simple  to  implement,  it  is  important  to  keep 
in  mind  that  it  is  an  indirect  method.  In  particular  this  procedure  (see  Fig.  2.2) 
minimizes  the  squared  magnitude  of  the  error 


eA[n]  —  x[n]  *  a[n]  —  b[n] 

where  the  sequences  a[n]  and  b[n]  are  defined  as 

def  J    an;    0  <  n  <  P    (a0  =  1) 


(2.13) 


a[n] 


0;      otherwise 


and 


,r  I  def  {  an\    0  <n  <  Q    (a0  =  1) 
0;       otherwise 


(2.14) 


(2.15) 


where  P  and  Q  are  the  number  of  poles  and  zeros  respectively.  Equation  2.13  repre- 
sents a  different  least  squares  error  from  that  in  (2.3)  where  the  quantity  to  minimize 
was  the  squared  magnitude  of  the  error  e[n]  (see  Fig.  2.1).  The  practical  significance 
of  this  difference  is  that  frequently  there  is  a  loss  of  accuracy  in  estimating  the  poles 
and  zeros  by  Prony's  method  [Ref.  3]. 

2.    Signal  Domain  Form  of  Prony's  Method 

An  alternative  formulation  of  the  method  described  above  can  be  obtained 
if  the  problem  (2.1)  is  stated  in  the  signal  domain  by  representing  x[n]  in  terms  of  a 
set  of  complex  exponentials: 


x[n]  «  C^"  +  C27"2  + h  CPTp 


(2.16; 


where  as  stated  before,  the  r*.  are  the  roots  of  the  polynomial  A{z),  assumed  to  be 
distinct,  and  the  complex  coefficients  Ck  provide  for  the  linear  combination  of  the  P 
roots. 

This  approximation  can  be  initially  formulated  as  in  the  section  above  and 
(2.9)  can  be  solved  in  the  least  squares  sense  for  the  coefficients  of  A[z)  (ie.  the 
vector  a).  The  roots  r*  of  A(z)  can  then  be  found  and  (2.16)  can  be  evaluated  for 
n  =  0, 1, . . . ,  iV,  —  1  to  produce  the  set  of  equations 
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(2.17) 


This  set  of  equations  can  then  be  solved  in  a  least  squares  sense  to  obtain  the  vector 
of  coefficients  c. 

In  the  case  of  multiple  roots  at  the  same  location,  a  slight  variation  of  the 
same  procedure  can  be  used.  Suppose,  for  example,  that  r\  is  a  double  root.  In  this 
case,  the  approximation  is 


x[n]  «  cxr*  +  c2nr\  + h  cPrP 


(2.18) 


and  the  matrix  equation  to  solve  for  the  coefficients  becomes 
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(2.19) 


This  situation  is  rare,  however,  because  computational  errors  and  errors  inherent  to 
the  modeling  method  itself  contribute  to  produce  roots  that  may  be  very  close  to 
each  other,  but  not  exactly  at  the  same  location. 


3.    Iterative  Prefiltering 

The  iterative  prefiltering  method  attempts  to  solve  the  ''direct"  problem 
mentioned  in  subsection  1  of  this  chapter  and  to  refine  the  initial  pole-zero  estimate 
by  solving  a  succession  of  linear  problems.  Equation  2.2  (error  for  the  direct  problem) 
can  be  written  in  the  ^-domain  as 

sW = *w  -  m = a-(2)  -  m = ^m^m.     (2,0) 

A[z)  A{z) 

The  notion  of  iteration  can  be  introduced  that  allows  computation  of  a  new  set  of 
poles  and  zeros  based  on  the  last  known  set  of  poles.  Iterative  prefiltering  replaces 
the  error  for  the  direct  problem  at  the  (i  +  l)th  iteration  with  the  iterative  error 

function 

/..,i            X(?)A{t+1)(z)  -  B(l+1)(z) 
g"+"(z)  =  A'-)/l      ^(z)" 52.  (2.21) 

If  /i^[rc]  is  the  inverse  ^-transform  of  1/A^(z),  then  the  error  can  be  written  in  the 
signal  domain  as 

e(,'+1)[n]  =  x[n]  *  h®[n]  *  a{l+1)[n]  -  b^l+1)[n]  *  h®[n].  (2.22) 

The  coefficients  a^!+1'[n]  and  6^+1^[n]  are  selected  to  minimize  the  corresponding  sum 

of  squared  errors 

Ns-i  , 

S(i+i)  =    £   \e(l+l)[n][  (2.23) 

n  =  P 

at  each  iteration;  this  situation  is  shown  in  Fig.  2.3.  No  general  proof  of  convergence 
has  been  given  for  this  algorithm;  however,  it  is  easy  to  see  that  if  the  iteration  does 
converge,  it  must  produce  the  same  answer  as  the  direct  method.  Specifically,  at 
convergence  A^  =  A^+l\  and  (2.21)  becomes  the  same  as  (2.20). 

If  we  use  an  indirect  modeling  procedure  like  Prony's  method  to  compute 
the  initial  vector  a  and  we  define  x^  as 

x^[n}  =  x[n]*h(t)[n],  (2.24) 
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Figure  2.3:  Block  diagram  for  the  iterative  prefiltering  method. 

then  the  sequences  /i^[n]  and  x^[n]  for  n  =  0, 1, . . . ,  Ns  —  1  can  be  computed  from 
the  recursive  difference  equations 

P 


hW[n]  =  6[n}-Y,<ik)h{l)[n-k} 


k=\ 


and 


x{l)[n]  =  x[n}-Yja^x{l)[n-k}. 


'2.25! 


(2.26) 


/t=i 


Thus  the  error  (2.22)  can  be  written  as 


e(,+1,M  =  £  4'+1V"[n  -  k]  -  £  if  l)k«[B  -  >]. 

j=0 


(2.27) 


fc=0 


In  order  to  find  the  coefficients  a£       and  bj       ,  the  error  is  written  in  matrix  form 


for  n  =  0,1,  ...,N,  -  1 
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and 
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.(t'+l) 


dp 


(•'+1) 


,       b^+1)  = 


6(f+l) 

°0 


l(*+1) 


(2.31) 


Equation  2.28  is  analogous  to  (2.10).    Thus  the  least  squares  problem  defined  by 
minimizing  the  norm  of  the  error  in  (2.28)  can  be  reduced  to 


([  x«   hw  ]*t[  xw  h«  ]) 


-  ^(i+i)  - 

a(.-+l)  ■ 

0 

b(i+l) 

.  0 

w 


here 


£(t+l)  _   ||e(t+l)|j2_ 

This  linear  equation  can  then  be  solved  for  the  vectors  of  filter  parameters. 


(2.32) 


(2.33) 


At  this  point,  it  is  instructive  to  compare  the  performance  of  the  three 
modeling  methods  outlined  in  this  chapter  by  applying  each  to  model  a  small  segment 
of  a  transient  sound  corresponding  to  a  wrench  being  struck.  This  wrench  sound  was 
recorded  and  sampled  in  the  laboratory  at  a  sampling  rate  of  approximately  10,240 
Hz.  This  signal,  denoted  by  wrenOl,  was  also  used  and  modeled  in  [Ref.  9].  Figure 
2.4  shows  a  100  point  segment  of  the  signal  wrenOl  with  one  of  the  three  models 
(Prony's,  signal  domain  form  of  Prony's,  and  Iterative  Prefiltering  )  overlaid.  In  all 
three  cases  the  signal  was  modeled  with  four  poles  and  four  zeros.  As  can  be  seen, 
the  difference  between  iterative  prefiltering  and  the  first  two  models  is  significant. 
The  non-iterative  methods  match  only  the  initial  points  of  the  sequence,  and  produce 
poor  approximations  in  modeling  the  remaining  part  of  the  signal.  This  is  due,  in 
large  part,  to  poles  not  sufficiently  close  to  the  unit  circle.  Iterative  prefiltering, 
on  the  other  hand,  produces  a  model  which  is  close  to  the  real  sequence  along  the 
complete  segment. 
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III.    ITERATIVE  ALGORITHM  IN  THE 
SIGNAL  DOMAIN 

A.  MULTIDIMENSIONAL  OPTIMIZATION  BY  GRADI- 
ENT METHODS 

A  multidimensional  function  f{x\,  X2,  ■  ■  •  ,  xn)  that  is  continuous  and  differen- 
tiable  can  be  minimized  using  one  of  several  very  powerful  hillclimbing  techniques 
known  as  gradient  methods  [Ref.  10,  p.  84].  Some  of  those  methods  are  derived  on 
the  basis  of  a  quadratic  model  that  can  be  obtained  from  a  truncated  Taylor  series 
expansion  of  /(x).  Let  x\k>  denote  the  value  of  x  at  the  A:th  iteration.  Then  for  any 
point  x  =  ySk'  +  6 ;  when  6  is  small,  the  function  can  be  approximated  by 

/(X«  +  S)  «  qW(6)  =  fW  +  g{k)T6  +  I*TG««  (3.1) 

where  g  and  G  represent  the  vector  of  first  derivatives  and  the  matrix  of  second 
derivatives  of  the  function  /(x)  respectively  and  they  should  be  available  at  every 
point.  In  Newton's  method  the  iterate  x^+1^  is  taken  to  be  x^fc'  -f  <5^  ,  where  the 
correction  6^  '  minimizes  q^(6).  This  method  is  only  well  defined  when  the  matrix 
of  second  derivatives  G  is  positive  definite,  in  which  case  the  Arth  iteration  of  Newton's 
method  is  given  by  the  following  procedure  [Ref.  11,  pp  44-46]: 

1.  solve  G^6  =  -gW     for    6  =  6(k) 

2.  setx(*+1)=x(fc)+5(A). 

(3.2) 

The  fact  that  G^  may  not  be  positive  definite  when  x^  is  far  from  the  solution, 
and  that  even  when  G^  is  positive  definite  convergence  may  not  occur,  makes  this 
method  undesirable  as  a  general  formulation  of  a  minimization  algorithm.  However,  a 
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number  of  variations  to  the  basic  method  have  been  proposed  that  are  more  suitable 
for  a  general  class  of  problems.  One  of  these  methods  is  Newton  s  method  with  line 
search  [Ref.  11,  pp.  47-49]  in  which  the  Newton  algorithm  is  used  to  generate  a 
direction  of  search 

s(*)  =  -GW-Vfc)  (3-3) 

which  can  later  be  used  in  a  line  search  algorithm  to  actually  calculate  the  correction 
6.  In  the  cases  when  G^k'  is  not  positive  definite,  the  linear  search  can  be  made 
along  ±s^  choosing  the  correct  sign  to  ensure  a  descent  direction.  However,  some 
difficulties  that  arise  here  (like  very  high  numerical  costs  and  failure  of  convergence 
for  some  special  cases)  make  this  an  undesirable  approach  for  our  algorithm. 

As  stated  before,  Newton?s  method  is  defined  only  when  the  matrix  G^  is 
positive  definite,  and  this  matrix  is  positive  definite  only  when  the  error  6  is  "small"; 
or  better  stated,  the  method  is  defined  only  in  some  neighborhood  tt(k>  of  x^  in 
which  q(k'(S)  agrees  with  /(x^  +  6)  in  some  sense.  In  such  cases,  it  is  correct  to 
choose  x(fc+1)  —  x(/c)  +  6(k\  with  the  correction  <5(/c)  minimizing  q^{6)  for  all  x(i)  +  6 
in  Q(k\  This  method  is  referred  to  as  the  restricted  step  method  because  the  step  is 
restricted  by  the  region  of  validity  of  the  Taylor  series.  [Ref.  11] 

The  region  of  definition  for  the  kth  iteration  can  be  expressed  as 

Ow  =  {x:  ||x-x(fc)||    <h(/c)}  (3.4) 

where  ||  •  ||  denotes  the  norm  of  the  vector.  In  this  case,  the  optimization  problem 
can  be  stated  as: 

minimize^  q{k){6)  subject  to  \\6\\    <  h{k) .  (3.5) 

As  mentioned  before,  the  least  squares  norm  is  the  one  most  commonly  used  in  this 
type  of  problems,  so  it  is  the  one  used  in  this  thesis  and  is  denoted  as  ||  •  ||2.  The 
problem  that  now  becomes  apparent  is  how  to  select  the  error  margin  h^)  of  the 


15 


neighborhood  (3.4).  This  margin  should  be  as  large  as  possible  because  the  iteration 
step  is  directly  related  to  it.  Various  methods  have  been  proposed  to  control  the 
parameter  h^;  one  of  these  methods  attempts  to  insure  that  the  Newton's  search 
direction  problem  (3.3)  is  always  defined  [Ref.  11.  pp.  100-103].  It  does  so  by  adding 
a  multiple  of  the  unit  matrix  I  to  G^'  and  computing  the  new  problem 

(G^  +  vl)6w  =  -g^  (3.6) 

where  the  net  effect  is  that  increases  in  v  cause  ||<$||2  to  decrease,  and  vice  versa. 
If  we  define 


r(*)  = 


A/(*)       /(*)_/(x(*>+*W) 


(3.7) 


AqW  fW-qW(6lk)) 

then  the  ratio  r^  represents  a  measure  of  the  accuracy  to  which  q^k\s^  ')  approxi- 
mates f(x^  +  <V  ')  on  the  A;th  step,  and  as  the  accuracy  increases  Ak^  gets  closer  to 
unity.  Using  (3.7),  Marquardt  [Ref.  12]  suggests  an  algorithm  that  tries  to  adaptively 
maintain  h^  as  large  as  possible  while  controlling  the  ratio  Ak\  The  kth  iteration 
of  such  an  algorithm  is  stated  as: 

1.  given  x(/c)  and  v(k\  calculate  g(/c)  and  G(fc); 

2.  factor  G^  +  i/"l;  if  not  positive  definite,  reset  v^  =  4i>^  and  repeat; 

3.  solve  (3.6)  to  find  6{k); 

4.  evaluate  /(x(/c)  +  6(k))  and  hence  r(/c); 

5.  if  rW    <   0.25  set  u^k+l)  =  4^k) 

else  if  r{k)    >   0.75  set  ^k+l)  =  v{k) /2 
else  set  iy(A:+1)  =  i/W; 

6.  if  rW    <   0  set  x^+1^  =  x<*>  else  set  x^+1)  =  x<*>  +  6(k). 

Here  the  parameters  0.25,  0.75,  4  and  2  are  arbitrary,  and  i/1)  >  0  is  also  chosen 
arbitrarily  [Ref.  11,  pp.  102-103].  Proofs  of  global  and  second  order  convergence  for 
this  algorithm  are  given  in  [Ref.  11,  pp.  96-98]  for  the  cases  when  the  first  and  sec- 
ond derivatives  of  the  function  f(x)  exist,  and  the  vector  x^)  belongs  to  a  bounded 
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n-dimensional  space  for  all  k.  Although  this  method  does  have  some  disadvantages, 
it  represents  a  good  basis  for  the  formulation  of  a  general  minimization  algorithm. 
Some  variations  of  this  method  were  considered,  but  it  was  found  that  for  the  spe- 
cific application  of  ARM  A  modeling  in  the  time  domain,  these  variations  produced 
an  extremely  high  overhead  in  calculations. 

B.  THE  ITERATIVE  PRONY  METHOD 

Let  us  now  return  to  the  problem  of  representing  a  sequence  x[n]  as  a  linear 
combination  of  complex  exponentials.  Equation  2.17  can  be  written  as 


Re  =  x  + 


(3.8) 


where  e  is  called  the  equation  error,  x  represents  the  data  which  may  or  may  not 
be  complex,  c  is  the  vector  of  complex  coefficients,  and  R  is  the  matrix  of  complex 
roots,  which  can  be  written  more  specifically  as 


R  = 


1  1 

rflj     +    jrh  rR2     +    jrh 

{rRl     +    jrhf  {rR2     +    jrhf 


1 
rRP     +    jrIp 
{rRP     +    jrIp): 


{rRp     +    jrIf 


\N.-i 


(3.9) 


„  (rRl     +    jrh)N>-1     (rR2     +    jn2)N^ 

where  rRi  and  rjt  represent  the  real  and  imaginary  components  of  the  i      root,  re- 
spectively, and  A^s  is  the  number  of  data  samples. 
By  defining 


def 


*T 


Qa=2  ||€||2  =  e*J6  =  (Rc-x)*i  (Rc-x), 


(3.10) 


it  is  clear  that  the  problem  is  to  find  the  vector  r  =  rR  +  j'17  of  P  complex  roots  that 
minimizes  the  function  Q(r). 

The  first  and  second  derivatives  of  Q  with  respect  to  r  are  represented  by  the 
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vector  g  and  matrix  G. which  are  defined  by 


and 


drIpdTRl  dTrpdrR2 


dQ 
drR2 


drRp 

do 
dQ 


dQ 


dri 


P     -i 


32Q 


drRl  drRp 
d2Q 


<32Q 


drR   3tj 


d2Q 


d2Q 


drRpdTRp 

d2Q 
dri18rRp 

d2Q 
dri2drRp 


82Q 


d2o 


9rR  ,3r/, 

32Q    ' 


dTR2dTRp        drR20rj1         drfydr^ 


d2Q 


TRpdrIl        drR    dri2 
d*Q  d*Q 


.'2Q 

dru  drri 

d2Q 


drj ,3ri2 

d2Q 
3r/23rjj         9rl2drl2 


d2Q 


dripdrRp 


d2Q 


92Q 


drjpdr^         drIpdri2 


(3.11) 


d2o 


drR  drip 

d2Q 


d2Q 
drRpdr,p 

d2Q 
dr[1dr[p 

32Q 
3r/.arjp 


d2Q 
dripdrif 


(3-12) 

In  order  to  provide  for  a  more  compact  notation,  define  the  gradient  operator  with 
respect  to  the  complex  vector  r  as  the  vector  of  partial  derivatives 


Vrd^f 


Vrfl 
Vr, 


8t 


«: 


9 


d 

3rj, 


drj. 


(3.13) 
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consequently,  (3.11)  can  be  expressed  as 


Equation  3.12  can  also  be  written  as 

dQ       do 

drRl    [    drRl        drp^ 


G  = 


dr 


«2 


92 


dQ 


drRl        drR2 


drRp 

d 

dri1 

d 


dQ 


dQ 


drRl        drR2 


■■<Q 

drRl 

dQ 

dTZ  [  drR\ 


dQ 


dQ 

dTR,, 

dQ 
drR2 


dQ 


dQ 

drRp 

dQ 
drRp 


dQ 
drRp 


dQ 


dr 


drIp    [    9tRi        drp2 


dQ 
drRp 


Vr*Q 

Vr/Q 


OQ 
dQ 


dri2 

dQ 

dri 


dQ         dQ 
drr,         drj. 


dQ 


'Hp  "'  l1 

dQ         dQ 
RP        drTl 


dQ 
drj2 

dQ 
drJr, 


dQ         dQ 
drIx        dTj2 


dQ 

drip 

dQ 

drj 


dQ 
dr,p 


dQ 
drip 

dQ 
drip 


dQ 
drip 


(3.14; 


(3.15) 


Now  using  a  somewhat  more  convenient  notation,  (3.15)  can  be  rewritten  as 


G 


ik  l 


dQ         dQ 


?rRt 


i  >r 


h 


d 
dri. 


\     dQ         dQ 
drRt        dr!t 


I  =  l...P 
k=l...P 


(3.16) 


and  it  becomes  clear  that  the  matrix  of  second  derivatives  can  be  expressed  compactly 


as 


G  =  VT 


(VrQ) 


71 


=  v, 


(VrfiQ)T    (Vr/Q)r 


(3.17) 


The  following  subsections  derive  explicit  expressions  for  the  two  quantities  g  and  G. 
1.    Vector  g  of  first  derivatives 


From  (3.8)  it  is  easy  to  see  that 
de        OR 


Or 


=  — — c    and 
dri 


de*T 
dr{ 


=  c 


dr{ 


(3.18) 
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where  r,  represents  the  real  or  imaginary  parts  of  the  ith  root.  Using  (3.18)  and  the 
chain  rule  in  (3.10)  leads  to 


3Q         mT  OR  xTdK'T 

C  +  C      — e. 


drRi  3rRt  drj^ 

and  explicit  evaluation  of  the  partial  derivatives  of  the  matrix  R  results  in 


(3.19) 


OQ 
9rRi 


=     e 


0 

1 

2(rRj+jr7i) 
3(rHj  +  >/,)' 


L  W-lJ^+jr/J 


AT. -2 


C8     + 


<f     0     1     2(rJJ|-<7>/<)    3(r* -jrz.)'    •••    (JV.  - 1)  (rj*  -jr,.) 


IV. -2 


(3.20) 


where  ct-  is  the  z     component  of  the  vector  c.  If  the  quantity  £,-  is  defined  as 


then 


def 


2(rfl,+jr/,)2 
3(rfli  +j>/,)2 


A',- 2 


Ct, 


_dR 


c    =    &, 


(3.21) 


(3.22) 


•  dR 


*r 


drfli 


and  (3.20)  can  be  written  compactly  as 

dQ 


drRt 


=  e      £i  +£;    e. 


(3.23) 


(3.24) 
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At  this  point  it  can  be  recognized  that  (3.24)  represents  the  addition  of  two  scalar 
quantities  where  the  first  one  is  the  complex  conjugate  of  the  second,  so  (3.24)  can 
be  further  simplified  to 


drRt 


=  2Rekfe 


(3.25) 


where  Re[-]  denotes  the  real  part  of  the  vector.   Finally,  using  (3.25)  for  i  =  1  ■  ■  ■  P 
in  the  upper  partition  of  (3.14)  produces 


VrfiQ  =  2Re< 


*r  -i 


T 


L  £p    J 


(3.26) 


A  similar  procedure  can  be  used  to  obtain  the  gradient  of  Q  with  respect 
to  the  imaginary  part  of  the  vector  r.  Once  more,  from  (3.18)  and  the  chain  rule 
applied  to  (3.10),  it  follows  that 

OR 


drIt 


,dR*2 
drh 


Jti 


-HI 


*T 


(3.27) 


(3.28) 


These  results  can  be  used  to  generate  an  expression  similar  to  that  in  (3.24)  for  the 
vector  of  first  derivatives  of  Q  with  respect  to  the  imaginary  components  of  r 

dQ 


drIt 


=  Jt*T^-jtf* 


(3.29) 


which  in  turn  defines  the  vector  of  partial  derivatives  with  respect  to  the  imaginary 
components  as 


Vr/Q  =  2Im<! 


if 

€     > 

{ [  if  \ 

j 

(3.30) 
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Equations  3.26  and  3.30  can  now  be  combined  as  shown  in  (3.14)  to  obtain  the  final 
expression  for  the  vector  of  first  derivatives 


VrRQ 

Vr/Q 


-  9 


Re<^ 


Ci 


dp  j 


Im  < 


f 

tf 

e    > 

Uf . 

j 

(3.31) 


2.   Matrix  G  of  second  derivatives 

An  expression  for  the  matrix  G  of  second  derivatives  can  be  obtained  as 
follows.  Substituting  (3.24)  and  (3.29)  into  (3.16)  yields 


G    = 


d 
drRk 

([•' 

■T«, 

+ efe 

J  [<'T(, 

)" 

d 

([•' 

T(, 

+<f<] 

J  [t"Tt,  ■ 

). 

1  =  1. 

..p 

k=  1 

..p. 

(3.32) 


Then  using  the  chain  rule  and  both  expressions  in  (3.18)  leads  to 
G  = 


drR,  dm. 


J 


.«T^i     |    r*T8R«rjf.  ,*T  dR  dC 

drR,    ^  C        arRt  5i         ««      arRtC  9rR] 


St  "r  «,■     dr     c  "I"    <9r 


9r/fc  dr^ 


I  =  1  .  .  .  P;       k-l...P. 


.*T?L 


+  C 


.J9R 


OLc         r'T  ^R 


I3R„       9C, 


yr/. 


(3.33) 
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From  the  definition  of  £  in  (3.21)  it  is  seen  that 


d£i  def 

S«Jfe    = 


drR 


6(rR,  +jrIt) 


Ns-l)(Na-2)(rRi+jrii 


Ns-3 


o, 


In  the  same  way  it  can  be  shown  that 

dif 


drRk 


sik 


C{,    if  i=k 


otherwise. 


(3.34) 


(3.35) 


drIk 


=    jslk 


(3.36) 


drIk 


=     -JSik 


(3.37) 


G 


These  last  four  expressions  together  with  (3.22),  (3. 23), (3. 27),  and  (3.28)  can  now  be 
substituted  in  (3.33)  to  obtain 

[e*TSlk  +  tfti  +  tftk  +  S^e]        J  [e*TStk  +  if  ^  -  tftk  -  age] 
J  [exTstk  -  tfti  +  tftk  -  age]     -  [e*Taik  -  ifti  -  l?U  +  sg"  e 

i  =  1 . . .  P 

k  =  l...P.  (3.38) 

Finally  notice  again  that  the  elements  of  the  matrix  G  are  formed  by  addi- 
tions and  subtractions  of  complex  scalars  with  their  respective  complex  conjugates; 
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thus  the  final  expression  for  G  is  given  by 


G    = 


2Rek:rd+2RekTel      21 


-21m 


t?tk   +2Im  sge      2Rckr^* 


Lm 


tftk 


-2Im|s*re 


-2  Re 


stk  e 


i  =  l. ..P;     j  =  l...P 


(3.39) 


where  it  should  be  remembered  that  s,-j.  =  0  for  all  i  ^  k. 
3.    Algorithm  implementation 

An  iterative  method  for  ARMA  modeling  in  the  time  domain  was  imple- 
mented using  the  results  of  the  last  two  subsections  in  conjunction  with  the  algo- 
rithm presented  in  section  A  of  this  chapter.  We  call  this  method  the  iterative  Prony 
method. 

The  algorithm  uses  the  signal  domain  form  of  Prony's  method  to  calculate 
an  initial  model  for  the  given  sequence.  From  there  it  uses  the  calculated  model 
to  compute  the  error  e,  the  vector  of  first  derivatives  g,  and  the  matrix  of  second 
derivatives  G  and  iterates  until  specific  conditions  are  met.  Figure  3.1  is  an  example 
of  how  the  algorithm  changes  the  position  of  the  poles  and  zeros  of  the  initial  model 
in  order  to  minimize  the  error.  This  figure  represents  the  poles  and  zeros  of  a  transfer 
function  of  order  (4,3) — 4  poles  3  zeros — that  was  overmodeled  using  a  (6,5)  order 
model.  It  is  clear  from  the  figure  that  the  tendency  in  this  case  is  to  have  a  pole-zero 
cancellation  (see  second  and  third  quadrants)  of  two  poles  and  two  zeros  as  expected. 

Some  features  were  added  to  the  basic  algorithm  in  order  to  deal  with  special 
modeling  cases.  Specifically,  if  the  initial  model  has  some  roots  on  the  real  axis,  then 
because  of  the  way  the  algorithm  iterates,  those  roots  never  move  away  from  the  real 
axis.  A  modification  was  therefore  introduced  to  deliberately  displace  those  roots 
from  the  real  axis  and  proceed  with  the  iterations.    If  the  tendency  of  the  roots  is 
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Figure  3.1:  Displacement  of  the  poles  and  zeros  of  an  iterative  Prony's  model 

"to  go  back"  to  the  real  axis,  then  they  are  returned  to  their  initial  position,  and  the 
iterations  continue.  Otherwise  the  roots  may  continue  to  spread  apart  and  move  as 
a  complex  pair.  Figure  3.2  is  an  example  of  the  displacement  of  the  poles  and  zeros 
of  an  order  (4,3)  model.  In  this  case  the  modeled  signal  actually  has  two  poles  on 
the  real  axis,  and  the  initial  model  correctly  placed  two  of  the  poles  in  the  real  axis. 
Those  poles  are  displaced  from  the  real  axis  by  the  algorithm,  but  then  after  some 
iterations  it  is  clear  that  the  poles  are  tending  to  return  to  the  real  axis.  At  this  point 
the  poles  are  forced  back  to  the  real  axis  by  setting  their  imaginary  parts  to  zero 
and  the  iterations  continue  until  convergence  is  obtained.  The  opposite  situation  is 
shown  in  Figure  3.3.  In  this  case  the  initial  model  also  has  two  poles  located  on  the 
real  axis,  but  contrary  to  the  case  presented  above,  the  roots,  after  being  displaced 
from  the  real  axis,  continue  to  move  away  from  the  axis  until  they  reach  their  final 
position  in  the  first  and  fourth  quadrants  closer  to  the  unit  circle. 
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Figure  3.2:  Displacement  of  the  poles  and  zeros  of  a  4-3  model.  The  modeled  signal 
in  this  case  actually  has  two  poles  on  the  real  axis 
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Figure  3.3:   Displacement  of  the  poles  and  zeros  of  a  4-3  model.   The  initial  model 
shows  two  poles  in  the  real  axis,  the  final  model  in  this  case  has  no  poles  in  the  axis 
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IV.    PERFORMANCE  OF  THE  ALGORITHM 
AND  MODELING  RESULTS 

A.  TEST  DATA  USED  IN  THIS  THESIS 

Two  types  of  signals  were  used  in  this  thesis  to  test  the  performance  of  the 
algorithm.  The  first  type,  which  will  be  called  simulated  test  data,  consists  of  five 
sequences  (tOl  to  t05)  each  one  hundred  points  long  that  were  produced  as  the 
impulse  response  of  a  known  rational  system.  Noise  to  produce  an  SNR  in  the  range 
of  10  to  15  dB  was  added  to  the  original  sequences,  and  the  resulting  sequences  were 
designated  as  t01_n  to  t05_n.  The  original  signals  are  described  in  Table  4.1  by  their 
transfer  functions  and  the  location  of  their  poles  and  zeros. 

The  second  group  of  test  signals  consists  of  recorded  acoustic  data.  Two  of  these 
signals  were  recorded  and  sampled  in  the  laboratory.  One  of  them  is  the  sequence 
wrenOl  already  mentioned  in  Chapter  II;  the  other  one  was  obtained  from  human 
speech,  in  particular,  the  signal  voweLa  corresponds  to  100  samples  of  the  Spanish 
vowel  a.  The  remaining  three  signals  from  the  group  of  acoustic  data  were  recorded 
at  sea  by  a  submarine  platform;  they  correspond  to  sounds  produced  by  marine  life 
and  ice  cracking.  The  description  of  the  acoustic  signals  is  presented  in  Table  4.2. 

B.  PRESENTATION  OF  RESULTS 

All  simulated  test  signals  were  modeled  twice  with  the  iterative  Prony  method. 
In  the  first  test  the  exact  number  of  poles  and  zeros  of  the  original  model  was  used; 
in  the  second  test  all  signals  were  modeled  using  two  more  poles  and  zeros  than  the 
original  model.  This  last  test  is  considered  closer  to  a  real  life  situation  where  the 
exact  order  of  the  signal  to  be  modeled  is  unknown. 
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TABLE  4.1:  DESCRIPTION  OF  SIMULA  TED  TEST  DATA 


NAME 


TRANSFER  FUNCTION 


POLES 


ZEROS 


tOl 


H(z] 


1-0.772-1 


l-i^gz^+o^osz-2 


0.9513  L±  0.6286 


0:  0.770 


t02 


H(z)  = 


0.744Z"1  -1.6993z-2+0.6302~3 
l-2.977z-1+3.909z-2-2.520z-3+0.717z-4 


0.9422  L  ±  0.6200 
0.8987  L±  0.6385 


0;  0.4686 
1.8069;  oo 


t03 


H(. 


0.171z-1-0.155z"3 


l+1.861z-1+2.588z-2  +  1.684z-3+0.819z-4 


0.9512  L  ±  1.8846 
0.9514  Z±  2.3041 


0;  0.9521 
-0.9521:  oo 


t04 


H(z 


0.981z~1-1.279z~2+0-871;~3 

l-0.956z-1+0.040z-2-0.705z-3+0.741z-4 


0.9513  L±  0.2097 
0.9049  L  ±  2.0943 


0:  oo 
0.9423  L  ±  0.8068 


t05 


H(?\   —      1-0.75Z"1 

v z/   l-i^Sz-'+o^gz-2 


0.8441;  0.9358 


0;  0.750 


TABLE  4.2:  DESCRIPTION  OF  ACOUSTIC  TEST  DATA 


SIGNAL 

DESCRIPTION 

wrenOl 

Transient  sound  corresponding  to  a  wrench  being  struck 

voweLa 

Spanish  vowel  a 

bio_2133a 

Sperm  whale 

bio_2385a 

Porpoise  whistle 

bio_80a 

Ice  cracking 
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To  mathematically  produce  a  meaningful  measure  of  the  performance  of  a  mod- 
eling algorithm  can  be  quite  difficult  since  various  norms  can  be  deceiving  when 
comparing  errors  of  signals  with  large  differences  in  magnitude.  Two  different  ap- 
proaches to  measure  the  performance  of  the  algorithms  were  therefore  used  in  this 
thesis.  The  first  is  quantitative  and  involves  computing  the  squared-norm  of  the  error 
between  the  model  and  the  actual  signal  and  dividing  it  by  the  total  energy  (norm) 
of  the  signal.  The  second  approach  is  to  overlay  in  a  plot  the  model  and  the  original 
signal  in  order  to  provide  a  visual  comparison  of  the  results.  This  is  less  quantitative 
but  frequently  more  revealing  of  errors  in  the  modeling  process. 

1.    Simulated  test  data 

The  first  data  sets  modeled  were  the  simulated  test  data  sets.  Figure  4.1(a) 
is  a  comparison  of  the  normalized  errors  that  result  when  the  sequence  t01_n  is 
modeled  with  2  poles  and  1  zero  using  both  iterative  prefiltering  and  iterative  Prony 
methods.  Figure  4.1(b)  and  (c)  show  100  points  of  the  sequence  t01_n  and  the 
two  order  (2,1)  models.  At  this  point  there  is  no  noticeable  difference  between  the 
iterative  prefiltering  and  the  iterative  Prony  models.  Figure  4.2(a)  again  shows  a 
comparison  of  the  normalized  errors  between  an  iterative  prefiltering  model  and  an 
iterative  Prony  model  of  t01_n  for  the  case  when  the  signal  t01_n  was  overmodeled 
using  models  of  order  (4,3).  Although  the  difference  between  the  two  modeled  signals 
in  this  case  is  not  large,  notice  that  the  error  for  the  iterative  prefiltering  method 
initially  increases  before  decreasing  while  the  error  for  the  iterative  Prony  method 
decreases  monotonically.  This  is  the  first  example  of  a  pattern  that  repeats  in  all  but 
one  of  the  simulated  test  signals  that  were  modeled.  Figures  4.3  through  4.10  give 
similar  comparisons  for  the  remaining  simulated  signals.  The  pattern,  which  can  be 
seen  in  Figures  4.1  to  4.10,  is  that  when  the  signals  are  modeled  with  a  number  of  poles 
and  zeros  different  from  that  of  the  actual  order  of  the  signal  (always  overmodeling 
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Figure  4.1:  Signal  t01.n  and  its  2  poles-1  zero  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t01-ii  and 
the  iterative  prefiltering  model,  (c)  Signal  tOLn  and  the  iterative  Prony  model. 
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Figure  4.2:  Signal  t01-n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t01-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t01-n  and  the  iterative  Prony  model. 
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Figure  4.3:  Signal  t02.n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t02-ii  and 
the  iterative  prefiltering  model,  (c)  Signal  t02-n  and  the  iterative  Prony  model. 
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Figure  4.4:  Signal  t02-n  and  its  6  poles-5  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  W2-n  and 
the  iterative  prefiltering  model,  (c)  Signal  W2-n  and  the  iterative  Prony  model. 
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Figure  4.5:  Signal  t03-n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  W3-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t.03-n  and  the  iterative  Prony  model. 
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Figure  4.6:  Signal  tOSja  and  its  6  poles-5  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t03-n  and 
the  iterative  prefiltering  model,  (c)  Signal  tOS.n  and  the  iterative  Prony  model. 
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Figure  4.7:  Signal  t04~n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t04~n  and 
the  iterative  prefiltering  model,  (c)  Signal  t04~n  and  the  iterative  Prony  model. 
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Figure  4.8:  Signal  t04~n  and  its  6  poles-5  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t04-n  and 
the  iterative  prefiltering  model,  (c)  Signal  W4-n  and  the  iterative  Prony  model. 
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Figure  4.9:  Signal  t05-n  and  its  2  poles-1  zero  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t05-n  and 
the  iterative  prefiltering  model,  (c)  Signal  t05-n  and  the  iterative  Prony  model. 
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Figure  4.10:  Signal  t05-n  and  its  4  poles-3  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  t05-n  and 
the  iterative  prefiltering  model,  (c)  Signal  W5-n  and  the  iterative  Prony  model. 
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in  this  case),  the  behavior  of  the  iterative  prefiltering  method  tends  to  degrade  (i.e., 
it  takes  longer  for  the  algorithm  to  reach  convergence)  while  the  behavior  of  the 
iterative  Prony  method  remains  the  same  or  even  improves  as  in  the  case  of  t05_n 
shown  in  Figures  4.9(a)  and  4.10(a).  Parts  (b)  and  (c)  of  Figures  4.1  through  4.10 
show  the  original  signals  and  their  respective  models  overlaid.  It  can  be  seen  that 
the  models  arrived  at  by  both  methods  follow  the  original  signals  very  closely  in  all 
cases.  Table  4.3  lists  the  location  of  the  poles  and  zeros  of  the  systems  used  to  model 
all  the  simulated  noisy  signals.  These  systems  were  obtained  using  both  the  iterative 
prefiltering  and  the  iterative  Prony  methods.  The  poles  and  zeros  shown  in  Table  4.3 
can  be  compared  to  the  poles  and  zeros  of  the  original  simulated  signals  presented  in 
Table  4.1.  It  is  clear  that  the  location  of  the  poles  and  zeros  of  the  modeled  signals 
should  not  be  exactly  the  same  as  those  of  the  original  signals  because  some  noise  (in 
the  order  of  10  to  15  dB  SNR)  was  intentionally  added  before  the  modeling  process. 
However,  Tables  4.1  and  4.3  show  a  close  relation  between  the  location  of  the  poles 
and  zeros  of  the  original  sequences  and  the  position  of  the  poles  and  zeros  of  the 
modeled  signals. 

2.    Acoustic  test  data 

As  mentioned  in  the  last  section  the  acoustic  test  data  represents  sounds 
recorded  both  underwater  and  in  a  laboratory  environment.  In  some  cases  shorter 
segments  were  selected  for  modeling  due  to  the  complexity  of  these  signals.  Once 
again  the  iterative  prefiltering  algorithm  was  used,  and  its  results  were  compared  to 
the  results  obtained  using  the  iterative  Prony  method. 

Before  presentation  of  results,  it  is  important  to  explain  how  the  model 
produced  by  the  iterative  prefiltering  method  was  selected.  For  all  the  cases,  the 
same  number  of  iterations  was  used  both  for  the  iterative  Prony  and  for  the  iterative 
prefiltering  methods.   In  order  to  obtain  the  best  possible  results  from  the  iterative 
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TABLE  4.3:   POLE-ZERO  LOCATION  OF  THE 
THE  SIMULATED  NOISY  TEST  DATA 


SYSTEMS  USED  TO  MODEL 


SIGNAL 

NAME& 
(ORDER) 

ITERATIVE 
PRONY 

ITERATIVE 
PREFILTERING 

POLES 

ZEROS 

POLES 

ZEROS 

t01_n 

(2,1) 

0.9511  Z±  0.6291 

0;  0.7612 

0.9507  L  ±  0.6290 

0;  0.7619 

t01_n 

(4,3) 

0.9510  Z±  0.6290 
0.9217  L±  3.1397 

0.7636 
0.9476;  0.9151 

0.9506  L±  0.6289 
0.9867  L  ±  2.0944 

0.7630 

0.9886  Z  ±  2.0880 

t02_n 

(4,3) 

0.9545  Z  ±  0.6238 
0.8455  L  ±  0.6743 

13.2102 
2.1568;  0.1799 

0.9250  L±  0.6313 
0.9142  Z±  0.6238 

6.8483 
1.6619;  0.4650 

t02_n 

(6,5) 

0.9541  Z  ±  0.6228 
0.8505  L  ±  0.6797 
0.9896  L  ±  2.7699 

17.1527 
2.1237;  0.0952 
0.8760  Z  ±  2.7876 

0.9301  Z±  0.6335 
0.9068  L±  0.6173 
0.9221  Z  ±  2.4366 

6.5260 

1.8203;  0.5542 
0.9800  Z  ±  2.3896 

t03_n 

(4,3) 

0.9510  Z±  1.8850 
0.9505  Z  ±  2.3029 

22.9950 
1.0511;  0.9365 

0.9499  Z±  1.8851 
0.9514  Z±  2.3030 

12.2468 
0.9680;  0.8020 

t03_n 

(6,5) 

0.9509  L  ±  1.8851 
0.9506  L  ±  2.3029 
0.9356  Z±  0.2779 

21.1093 
0.9340;  0.6882 
1.1445  L  ±  0.3507 

0.9492  Z  ±  1.8854 
0.9515  L±  2.3028 
0.9796  Z±  1.5986 

7.5747 

0.9566;  0.6099 
0.9870  Z  ±  1.6035 

t04_n 

(4,3) 

0.9022  Z  ±  2.0981 
0.9515  L±  0.2097 

22.9853 

0.9494  L±  0.8212 

0.9048  Z  ±  2.0956 
0.9516  Z±  0.2090 

2362.7 

0.9000  Z±  0.8090 

t04_n 
(6,5) 

0.9040  Z  ±  2.0966 
0.9517  Z±  0.2094 
0.7237;  0.7221 

17.5349 

0.9404  L  ±  0.8333 

0.7099  L±  3.1147 

0.9053  Z  ±  2.0933 
0.9518  L±  0.2093 
0.9542;  0.8105 

57.6799 

0.9486  Z  ±  0.8189 

0.9461;  0.7803 

t05_n 

(2,1) 

0.9416;  0.7840 

0.6693 

0.9363;  0.8325 

0.7236 

t05_n 

(4,3) 

0.9366;  0.8345 
0.9571  Z±  2.8523 

0.7330 

0.9571  L±  2.8436 

0.9368;  0.8231 
0.8839  Z  ±  2.7995 

0.7019 

0.8998  Z  ±  2.7977 
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prefiltering  algorithm  (especially  in  those  cases  when  convergence  was  not  obtained), 
the  model  selected  was  that  corresponding  to  the  iteration  with  the  lowest  error 
between  the  model  and  the  original  sequence.  Figure  4.11  is  an  example  of  one  of 
the  cases  when  convergence  was  not  reached  using  the  iterative  prefiltering  method. 
It  can  be  seen  that  if  the  system  selected  is  the  one  obtained  at  the  last  iteration 
(20t/l)  of  the  algorithm  (which  in  this  case  corresponds  to  a  peak  of  the  error),  then 
the  resulting  model  does  not  follows  the  original  signal  at  all.  However,  if  we  select 
the  system  from  the  iteration  for  which  the  error  is  the  smallest  (17t/l  iteration  in  this 
case),  then  we  obtain  a  model  that  is  closer  to  the  original  signal.  Some  insight  can 
also  be  obtained  if  we  look  at  the  behavior  of  the  poles  and  zeros  of  the  system  while 
the  error  of  the  iterative  prefiltering  algorithm  is  oscillating.  Figures  4.12  and  4.13 
show  the  position  of  the  poles  and  zeros  during  the  oscillation  that  exists  between 
iterations  15  to  20  of  the  case  presented  in  Figure  4.11.  For  the  first  part  of  these 
iterations  when  the  error  is  low,  the  poles  remain  at  approximately  the  same  position. 
At  iteration  20,  (Fig.  4.12(f))  the  poles  suddenly  spread  apart,  some  to  even  outside 
the  unit  circle,  causing  the  large  increase  in  the  error.  A  similar  effect  occurs  with 
the  zeros,  although  there  is  some  significant  movement  of  the  zeros  in  the  earlier 
iterations.  This  type  of  behavior  is  avoided  in  the  iterative  Prony  method  where  a 
controlled  and  systematic  displacement  of  the  poles  and  zeros  is  produced  to  finally 
reduce  the  error  between  the  modeled  and  original  signals. 

The  same  approach  used  in  the  last  subsection  was  used  to  model  the  acous- 
tic data.  All  signals  were  first  modeled  with  a  low  order  system  and  subsequently 
remodeled  using  a  higher  order  system.  The  results  are  shown  in  Figures  4.14  to  4.23. 
In  most  of  the  cases  the  models  obtained  using  the  iterative  Prony  method  closely 
follow  the  original  signals.  It  can  also  be  observed  that  the  iterative  Prony  method 
does  not  have  as  large  a  dependency  in  the  order  of  the  model  selected  as  does  the 
iterative  prefiltering  method.  While  it  is  of  course  necessary  to  select  a  model  with  a 
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Figure  4.11:  Signal  voweLa  being  modeled  by  iterative  prefiltering  using  a  (10,9) 
order  system,  (a)  Normalized  squared-norm  of  the  error  between  the  model  and  the 
actual  signal,  (b)  Signal  voweLa  and  the  iterative  prefiltering  model  selected  from 
the  20t/l  iteration,  (c)  Signal  voweLa  and  the  iterative  prefiltering  model  selected 
from  the  17t/l  iteration. 


43 


1 

\ 

X 

0 

X 

X 
X 

x* 

-1 

Re[] 

-2 

(a) 


(b) 


1 

1 

0 

X 

— x '•■■'  '■'         x  •          - 

X 

i— * 

0 

X 

X                  •                   X 
X 

■1 

Re[] 

-1 

Re[] 

2 

-2 

(c) 


(d) 


? 


Re[] 


lh 
0 
-1 

-2 


Xx 

X        X 

X 

X        X 


Re[] 


(e) 


(f) 


Figure  4.12:  Behavior  of  the  poles  of  the  system  used  to  model  the  signal  voweLa 
during  one  oscillation  of  the  iterative  prefiltering  algorithm,  (a)  15</l  iteration,  (b) 
16th  iteration,  (c)  17</l  iteration,  (d)  18t/l  iteration,  (e)  19t/l  iteration,  (f)  20£/l 
iteration. 
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Figure  4.13:  Behavior  of  the  zeros  of  the  system  used  to  model  the  signal  voweLa 
during  one  oscillation  of  the  iterative  prefiltering  algorithm,  (a)  15th  iteration,  (b) 
16th  iteration,  (c)  17i/l  iteration,  (d)  18f/l  iteration,  (e)  19i/l  iteration,  (f)  20t/l 
iteration. 
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high  enough  order  to  properly  model  the  original  signal,  once  such  a  model  has  been 
selected,  increments  or  small  variations  to  its  order  do  not  produce  degradations  on 
the  performance  of  the  algorithm.  On  the  contrary,  it  was  found  that  in  some  cases 
the  performance  of  the  iterative  prefiltering  algorithm  can  be  significantly  reduced 
when  the  order  of  the  model  is  increased  slightly.  It  can  also  be  seen  that  the  itera- 
tive Prony  algorithm  tends  to  provide  a  closer  match  to  the  data  than  the  iterative 
prefiltering  method,  and  in  most  of  the  cases  the  rate  of  convergence  of  the  iterative 
Prony  method  was  higher  than  that  of  iterative  prefiltering.  Another  important  point 
that  can  be  extracted  from  the  results  presented  in  Figures  4.14  to  4.23  is  that  while 
convergence  with  neither  algorithm  is  guaranteed,  we  obtained  convergence  with  the 
iterative  Prony  method  in  all  cases  for  these  acoustic  signals.  The  same  was  not 
true  for  iterative  prefiltering.  In  most  of  the  cases  the  error  for  the  iterative  Prony 
method  begins  to  decrease  starting  at  the  first  iteration,  and  although  the  change  is 
not  monotonic  in  all  cases,  the  error  after  a  few  iterations  is  consistently  lower  than 
the  initial  error. 
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Figure  4.14:  Signal  wrenOl  and  its  7  poles-6  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  wrenOl  and 
the  iterative  prefiltering  model,  (c)  Signal  wrenOl  and  the  iterative  Prony  model. 
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Figure  4.15:  Signal  wrenOl  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
wrenOl  and  the  iterative  prefiltering  model,  (c)  Signal  wrenOl  and  the  iterative 
Prony  model. 
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Figure  4.16:  Signal  voweLa  and  its  10  poles-9  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
voweLa  and  the  iterative  prefiltering  model,  (c)  Signal  voweLa  and  the  iterative 
Prony  model. 
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Figure  4.17:  Signal  voweLa  and  its  14  poles-13  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
voweLa  and  the  iterative  prefiltering  model,  (c)  Signal  voweLa  and  the  iterative 
Prony  model. 
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Figure  4.18:  Signal  bio.2133a  and  its  8  poles-7  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio.2133a  and  the  iterative  prefiltering  model,  (c)  Signal  bio.2133a  and  the  iterative 
Prony  model. 
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Figure  4.19:  Signal  bio.2133a  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio.2133a  and  the  iterative  prefiltering  model,  (c)  Signal  bio.2133a  and  the  iterative 
Prony  model. 
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Figure  4.20:  Signal  bioJ2385a  and  its  8  poles-7  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bioJ2385a  and  the  iterative  prefiltering  model,  (c)  Signal  bioJ2385a  and  the  iterative 
Prony  model. 


54 


— 

a 

— 

4> 


1U' 

—  r 

-J 

i                      i 

L— 

_ 

10° 

Z 

solid  =  iterative  Prony  method 

= 

dashed  =  iterative  prefiltering  method 

~- 

10-1 

1 

- 

m-2 

i                      i                     i 

Iteration  number 

- 

xlO4 


10 


15 
(a) 


20 


25 


30 


solid  =  signal  bio_23  85a 

dashed  =  iterative  prefiltering  model 


0 


50 


xlO4 


solid  =  signal  bio_2385a 
dashed  =  iterative  Prony  model 


Figure  4.21:  Signal  bioJ2385a  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bio.2385a  and  the  iterative  prefiltering  model,  (c)  Signal  bio.2385a  and  the  iterative 
Prony  model. 
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Figure  4.22:  Signal  bioSOa  and  its  8  poles-7  zeros  models,  (a)  Normalized  squared- 
norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal  bioSOa  and 
the  iterative  prefiltering  model,  (c)  Signal  bioSOa  and  the  iterative  Prony  model. 
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Figure  4.23:  Signal  bioSOa  and  its  12  poles-11  zeros  models,  (a)  Normalized 
squared-norm  of  the  error  between  the  models  and  the  actual  signal,  (b)  Signal 
bioSOa  and  the  iterative  prefiltering  model,  (c)  Signal  bioSOa  and  the  iterative 
Prony  model. 
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V.    CONCLUSIONS 

A.   DISCUSSION  OF  RESULTS 

In  this  thesis,  a  new  method  for  modeling  signals  in  the  time  domain  is  developed 
and  applied  to  model  both  recorded  acoustic  data  and  simulated  signals  produced 
as  the  impulse  response  of  a  known  system.  We  call  this  method  the  iterative  Prony 
method.  In  most  of  the  simulated  test  data  sets  the  models  provided  by  the  iterative 
Prony  method  are  sufficiently  close  to  the  original  signals;  in  most  cases,  it  is  diffi- 
cult to  distinguish  between  the  signal  and  the  model.  When  modeling  the  acoustic 
data  distortion  becomes  apparent  in  some  of  the  models,  which  may  be  due  to  the 
complexity  of  the  structure  of  the  signals.  However,  this  distortion  is  no  worse  than 
for  any  of  the  best  algorithms  that  have  been  used  to  model  this  data  previously. 

The  new  algorithm  was  compared  very  specifically  to  the  iterative  prefiltering 
algorithm  [Ref.  7,  8]  which  has  been  used  in  modeling  a  variety  of  acoustic  data 
[Ref.  3,  9].  The  rate  of  convergence  of  the  iterative  Prony  method  was  in  most  of  the 
cases  comparable  or  superior  to  that  of  the  iterative  prefiltering  algorithm.  Thus, 
while  iterative  prefiltering  sometimes  has  convergence  problems,  the  new  algorithm 
is  much  more  dependable  in  that  respect.  The  price  to  pay  for  this  improvement  is 
in  the  number  of  computations.  While  the  number  of  floating  point  operations  per 
iteration  in  the  iterative  prefiltering  method  is  approximately 

64(P  +  Q  -  l)3  +  8NS(P  +  Q)2  +  10(P  +  Q)NS  +  12NS, 

iterative  Prony  requires  about 

672P3  +  (24JV,  +  102)P'2  +  {60NS  +  46)P  +  198iVs 
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floating  point  operations  at  each  iteration.  For  example  for  a  complex  data  set  of 
length  100  (Ns  =  100)  and  P  =  Q  =  8  we  have  approximately  452,400  floating 
point  operations  per  iteration  using  iterative  prefiltering  versus  572,360  using  iter- 
ative Prony.  If  we  increase  the  order  of  the  model  to  P  =  Q  —  16,  then  we  have 
approximately  2,787,824  operations  per  iteration  in  the  iterative  prefiltering  algo- 
rithm versus  3,509,560  in  the  iterative  Prony  method. 

B.   RECOMMENDATIONS  FOR  FUTURE  WORK 

The  iterative  prefiltering  algorithm  has  been  the  main  tool  in  the  modeling 
efforts  for  sonar  data  modeling  [Ref.  9,  13].  The  new  iterative  Prony  algorithm  is 
now  at  a  stage  where  it  can  be  substituted  for  the  iterative  prefiltering  algorithm  and 
tested  in  operational  use.  To  do  so  needs  some  further  programming  to  make  the 
segmentation  of  the  data  automatic  and  to  make  the  entire  modeling  procedure  more 
of  a  "turn  crank"  operation.  These  should  be  some  of  the  very  next  steps.  In  addition, 
the  practical  implications  of  the  increased  computation  needs  to  be  addressed,  and 
if  possible  new  methods  need  to  be  developed  to  help  reduce  computations. 

In  a  larger  sense  the  work  reported  in  this  thesis  can  be  used  as  a  base  for 
possible  applications  of  the  iterative  Prony  method  in  the  problems  of  filter  design, 
speech  processing,  and  spectral  estimation.  The  expressions  for  the  vector  of  first 
derivatives  g  and  the  matrix  of  second  derivatives  G  of  the  error  derived  in  Chapter 
III  can  be  used  along  with  different  minimization  methods  to  provide  for  other  new 
modeling  methods  that  may  adapt  better  to  specific  modeling  problems. 
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